Qunatifying the effect of extrinsic factors in noisy gene expression 



Purushottam D. Dixit 
Biosciences Department, Brookhaven National Laboratory, Upton NY 1197^ 

We present a maximum entropy framework to separate the intrinsic and the extrinsic contribu- 
tions to noisy gene expression solely from the profile of expression. We express the experimentally 
accessible histogram of mRNA copy numbers by accounting for possible variations in global extrinsic 
factors. The distribution of extrinsic factors is estomated using the maximum entropy principle. 
Our results show that extrinsic factors quantitatively and qualitatively affect the histogram. Specifi- 
cally, we suggest that the cell to cell variation in extrinsic factors accounts for the wider that Poisson 
distribution of mRNA copy numbers. We successfully test our framework on a numerical simulation 
of a simple gene expression scheme that accounts for the variation in extrinsic factors. We also make 
falsifiable predictions, some of which are tested on previous experiments in E. coll while others need 
verification. Application of the current method to understanding noise in mRNA production in 
eukaryotes and protein production is also discussed. 
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Recent experiments show that the life cycle of a gene 
product inside the cell is stochastic. For any gene, there 
exists great cell to cell variation in the expression level 
of both the protein and the mRNA (jTHT0[) and changing 
this variation has phenotypical effects (fTTHT4|) . Recently, 
it was also shown that co-regulated proteins have corre- 
lated variability (|T5|) . This variation arises from a) the 
'intrinsic' statistical mechanical fluctuations in diffusion 
and binding of the molecules involved in gene expression 
and b) the variation in global 'extrinsic' factors such as 
levels of various molecules such as RNA polymerase, ri- 
bosome, proteases, and RNAses that differ from cell to 
cell dU[16]). I n a given population of cells, the noise 



(m 2 ) — (m) 



serves as a useful quantification of variability in gene ex- 
pression where (m) is the mean level of the gene product 
m (mRNA or protein) and (m 2 ) — (m) 2 is the variance. 

Experimentally, the contribution to r/^ associated with 
extrinsic factors, the extrinsic noise 7Je, can be measured 
separately from the intrinsic noise rji (0 [J5J HU) . It 
is now known that, for proteins, the extrinsic noise r/E is 
the dominant contributor to 7/t © and can change the 
profile of gene expression in a non-trivial manner (|17p . 
Evidently, an important step towards conceptual under- 
standing of the noisy gene expression is to quantitatively 
account for the extrinsic cell to cell variations in the ma- 
chinery of transcription and translation. To the best of 
our knowledge, there exists no comprehensive framework 
to do so. 

The major technical hurdle in building a comprehen- 
sive theory for extrinsic noise originates in the multitude 
of factors that contribute to it. Instead of accounting 
for all the extrinsic contributors ab initio, we develop 
a maximum entropy framework to estimate their effect 
on gene expression from a few experimental observables. 
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The framework allows us to quantify the cell to cell in 
the extrinsic factors. We successfully test our results on 
a simple numerical scheme that incorporates these ef- 
fects. Most importantly, we show that extrinsic factors 
can quantitatively and qualitatively affect the experimen- 
tally observed histogram of the gene expression product 
(protein or mRNA). 

In this work, we prefer an analytical treatment and fo- 
cus our attention to the less explored variation in mRNA 
copy numbers in prokaryotes. Generalization of the cur- 
rent framework to understanding noise in mRNA pro- 
duction in eukaryotes and protein production is discussed 
and left for further studies. 



I. METHODS 
A. The distribution of mRNA copy numbers 

Consider the simplified reaction scheme 

DNA -k mRNA A 4> 

of transcription and degradation of mRNA molecules of 
a particular gene. 7 is the rate of transcription and 5 is 
the rate of degradation. We have neglected the activa- 
tion states of the DNA molecule e.g. promoter fluctua- 
tions and chromatin remodeling (fH 151 [TO)): an excellent 
assumption in studying prokaryotic gene expression. The 
steady state solution of the reaction scheme is the Pois- 

— k m 

son distribution P(m\k) — e — f of mRNA copy number 
m ([P8[) with effective synthesis rate k — j/S. 

Notice that the effective synthesis rate k depends in 
a complicated manner on the states of many molecules 
in the cell including the components of RNA polymerase, 
the dynamics of assembly of the RNA polymerase holoen- 
zyme, various RNAse molecules, and other competing 
genes. Consequently, it varies from cell to cell and also as 
a function of time from the start of the cell cycle. Thus, 
while studying gene expression in a population, instead of 
fixing a particular value of the effective synthesis rate k, 
we need to consider P(k) the probability distribution of 
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k. P(k) quantifies the extrinsic contribution noisy gene 
expression. 

For a given gene, experimentally assessing the cell to 
cell variability in k is non-trivial — P(k) has to be in- 
ferred from limited experimental information viz. mean 
expression level, cell to cell variation in gene expression 
level, etc. Here, we use the maximum entropy formal- 
ism to estimate P(k) using only the mean gene expres- 
sion level. Using the method of Lagrange multpliers, we 
restrict the maximization of entropy to candidate func- 
tions P*{k) that reproduce the experimentally observed 
mean expression. See supplementary materials for de- 
tails of all derivations. We direct interested readers to 
ref. (|19j) which treats a generalized maximum entropy 
problem similar to the one dealt here. We find that the 
functional form of P(k) is a Gamma distribution, 



(1) 



In Eq. 16 /it is the mean expression level and a = ?7i/?7e 



is the ratio of the intrinsic and the extrinsic noise. The 
joint distribution P(m, k) is then given by, 

g— ak ym+iict— 1 



P(m, k) = P(m\k)P(k) cx 



(2) 



The experimentally accessible histogram P(m) is ob- 
tained by summing over all variations in k i.e. summing 
over the effect of cell to cell variation in extrinsic factors, 



P(m) = ^ P( m , k)oc^2 
k 



We estimate -P(m) to be the negative binomial distribu- 
tion (the discrete version of the gamma distribution) , We 
estimate -P(m) to be the negative binomial distribution 
(the discrete version of the gamma distribution), 



P(m) oc 



1 



r [m + afj] 



(1 + aY 



(3) 



B. Noise decomposition of experimental data 



We estimate the total noise rj^ from Eq. 17 (see supple 
mentary materials for details), 



and, 
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The greater than Poisson relationship between t]t and 
the mean mRNA copy number /i (see Eq. 20 1 is some- 
times attributed to non-Poissonion dynamics e.g. pro- 
moter fluctuations, chromatin remodeling, mRNA syn- 
thesis bursts etc. flU H OH. We suggest that this 
greater than Poisson relationship may be present even 
in the absence of such complicated dynamics and could 



be attributed to the cell to cell variation in other extrin- 
sic factors. In fact, experiments show that in E. coli (|7|) . 
a system where promoter fluctuations and chromatic re- 
modeling may not be important, the noise in mRNA lev- 
els is still greater than expected from a Poisson process 
and the distribution of mRNA copy numbers wider than 
expected from a Poisson process. 

If exotic DNA dynamics are known to be not im- 
portant, Eq. [17] and Eq. [20] completely determine the 
histogram of mRNA copy numbers from experimentally 
measured mean expression level fi and total noise i]t. 
Moreover, r/i and 77E can be estimated from the histogram 
alone. Importantly, the framework estimates the hitherto 
illusive effect of cell to cell variation in extrinsic factors 
on gene expression viz. the distribution P(k) of the ef- 
fective synthesis rate k. Moreover, the joint distribution 
(Eq. [2]) allows us to estimate joint moments for exam- 
ple the covariance ^^/{'fc) ^ ' which we predict to be 
equal to the extrinsic noise t/e (see supplementary mate- 
rials) . These are some of the falsifiable predictions of the 
development presented here. 



II. RESULTS AND DISCUSSIONS 

A. Numerical validation 

We analyze a simple numerical scheme of mRNA syn- 
thesis that accounts for variations in extrinsic factors us- 
ing the development presented here (see Fig. [I] for a car- 
toon supplementary materials for details). Briefly the 
rate of synthesis 

7 = 7o [RNAP] 

and the rate of degradation 

S = 6 [RNAase] 

of rGene, the mRNA of the gene under consideration, 
both depend on the concentration of the cellular proteins 
that carry out those reactions viz. [RNAP] (a proxy for 
the RNA polymerase complex) and [RNAase] (a proxy 
for RNAase) respectively. Both the proteins are them- 
selves are stochastically synthesized and degraded. The 
variation in the proxies mimics the cell to cell variations 
in extrinsic factors. The effective synthesis rate k is di- 
rectly proportional to the ratio [RNAP] /[RNAase]. We 
implement the Gillespie algorithm ([20]) to estimate the 
steady state distribution of mRNA copy numbers from 
the scheme. 

In order to clearly elucidate the effect of extrinsic fac- 
tors on gene expression profile, in Fig. [2] we show the his- 
togram of mRNA copy numbers for three different levels 
of noise (quantified by the coefficient of variation in 
k) keeping the mean expression constant. 

In the left panel of Fig. [2] we show the histogram of 
mRNA copy numbers when the coefficient of variation 
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FIG. 1. A cartoon of the simplified scheme of mRNA production that takes into account extrinsic factors in gene expression 
levels (see supplementary materials for details). In the scheme, RNAP serves as the proxy for the RNA polymerase holoenzyme 
complex, RNAase is the proxy for RNA degradation machinery. The rate of synthesis of rGene, the RNA of a given gene is 
directly proportional to the concentration [RNAP] of the protein product of the RNAP gene. Similarly, the rate of degradation 
of rGene is directly proportional to the concentration [RNAase], the protein product of RNAase gene. RNAP and RNAase 
themselves are synthesized and degraded stochastically. 
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FIG. 2. The histogram of mRNA copy numbers (red dots), the Poisson distribution fit (black dashes) and the marginal 
distribution fit (solid blue, see Eq. 171 for three different scenarios in the numerical simulation. The mean mRNA copy number 
fx ss 4.4 for all three cases. Left: small variations in extrinsic factors (r/k ~ 5 x 10 -5 ) results in a histogram of mRNA copy 
numbers that is well described by a Poisson distribution. Middle: higher variation in extrinsic factors (rjk ~ 2.5) broadens the 
histogram of mRNA copy numbers. The marginal distribution -P(m) (see Eq. |17[ ) fits the data well (compare with right panel 
of Fig. Q. Right: high variation in extrinsic factors (r]h ~ 3.8) . Again, note that the histogram of mRNA copy numbers is 
wider than a Poisson distribution and the marginal distribution P (m) fits the simulation well. 



T}k, defined as 



(fc 2 ) - (k)' 
(k) 2 



, is low (r)k ~ 5 x 10 -5 ). Observe that the histogram of 
mRNA copy numbers (red circles) is well described by 
a Poisson distribution (black dashes), as is expected. If 
we increase the variation in k (rjk ~ 2.5 in the middle 
pannel and rjk ~ 3.8 in the right panel), the histogram of 
mRNA copy numbers gets broader and is best described 
by P(m) (Eq. 17 solid blue) rather than Poisson distri- 
bution (dashed black). Thus, even though the mRNA 
synthesis and degradation is governed by a Poisson pro- 
cess with an effective synthesis rate fc, the cell to cell 
variation in the rate itself makes gene expression a dou- 
bly stochastic process (fTT)l [2"Tj) and leads to a histogram 
of mRNA copy numbers that is not Poisson-distributed 
and is best described by a Gamma-like distribution. 
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FIG. 3. The observed histogram of mRNA copy numbers for 
the gene TufA © is wider than a Poisson distribution (dashed 
black) . Our results predict that the experimentally measured 



mRNA copy number histogram is described by Eq. 17 (solid 
blue). r\k ~ 0.7 is the estimated coefficient of variation of the 
effective synthesis rate k. 



B. Interpreting experiments 

We analyze the experimental data on mRNA abun- 
dance from Taniguchi et al. ([7]) and interpret it within 



the current framework. 

Fig. [3] shows the best fit to the histogram of mRNA 
copy numbers for the E. coli gene TufA ([7]). The Poisson 
distribution does not capture the mRNA histogram while 
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FIG. 4. Left: Log binned average trends in the intrinsic (□), extrinsic (•), and total (a) noise (3, 6, 16) as a function of fi, 
the mean mRNA level. Total noise and mRNA copy numbers were taken from Taniguchi et al. (0 and the decomposition into 
intrinsic and extrinsic noise is estimated from Eq. |20| The intrinsic and the extrinsic noise both contribute significalntly in 
total noise. The extrinsic noise saturates at high expression levels (red circle) (fll l4l 171 l8l I15[). Right: The estimated coefficient 
of variation rjk of the effective synthesis rate k as a function of mean expression level p (see Eq. 16 1. The scatter plot and the 
running average (red line) suggest that there exists a lower limit to r)k- 



Eq. 17 describes it well (for a comparison with numerical 
simulations, see the right panel of Fig. [2]). See sec. Ill for 
a discussion about the goodness of fit. 

In the left panel of Fig. |4j we show the predicted log- 
binned average trends in the decomposition of the total 
noise into its intrinsic and extrinsic components. The 
components are estimated from Eq. [20] for ~ 130 genes 
reported in Taniguchi et al. 0. The noise decreases as 
mean expression level increases and both intrinsic and 
extrinsic components contribute significantly to the total 
noise. The total noise and the extrinsic noise saturate at 
high expression levels, sometimes referred to as the 'ex- 
trinsic limit' ([41 [71 151 IT5|). Importantly, our framework 
also allows us to directly estimate the cell to cell varia- 
tion P(k) of the effective synthesis rate k. In the right 
panel, we show that indeed as the average expression level 
of mRNA fj, increases, the coefficient of variation de- 
creases and appears to saturate at high expression levels 
clearly elucidating the extrinsic limit to the noise in the 
expression level of genes. 

The mRNA histogram from a slightly involved 
model that captures the activation state of the DNA 



molecule ([10]) results in a distribution identical to Eq. 17 
In Eukaryotic systems, where promoter fluctuations, 
chromatin remodeling, and mRNA synthesis bursts are 
known to be important (pH f5T fTU|) : one can generalize the 
current framework to obtain the distribution of mRNA 
copy numbers that includes the variation in extrinsic 
noise. Application of the current framework to the noise 
in protein production is also straightforward. Unfortu- 
nately, the distribution of copy numbers in the above two 
cases is described by two parameters. The development 
presented here adds one additional parameter to describe 
the extrinsic noise. Thus, the resulting distribution will 
be described by three parameters. Consequently, analyz- 
ing the currently reported experimental measurements of 
noise in mammalian mRNA levels and/or protein levels 



will be an overfit. We leave it for further studies. 



III. CONCLUDING REMARKS 

Measurements of cell to cell variation in protein num- 
bers show that the extrinsic contributions play a dom- 
inant role ([3]). Yet, much of the theoretical develop- 
ment in understanding noise in gene expression has fo- 
cussed on the effect of intrinsic contributors viz. sta- 
tistical mechanical fluctuations in binding and diffusion 
of molecules. The limited treatment extrinsic noise 
has recieved ([7J [151 employs the linear fluctuation- 
dissipation like susceptibility analysis (|16p or ad hoc as- 
sumptions about the nature of variation in extrinsic pa- 
rameters ([7] l2"2]) . 

To the best of our knowledge, we have, for the first 
time, presented a framework that systematically sepa- 
rates the intrinsic and the extrinsic contributors to noisy 
gene expression from limited information about the gene 
expression profile. We conclude that extrinsic factors can 
quantitatively and qualitatively change the experimen- 
tally accessible histogram of mRNA copy numbers. Im- 
portantly, the framework allows us to directly estimate 
the hitherto elusive variation in global extrinsic factors. 

Specifically, we show that even if mRNA synthesis and 
degradation is described by a simple Poisson process, ow- 
ing to the cell to cell variation in the effective synthesis 
rate k, the experimentally accessible histogram of mRNA 
copy numbers is broader and we estimate it to be the neg- 
ative binomial distribution (see Eq. 17). Consequently, 



we find that the greater than Poisson relationship be- 
tween noise 7/t and the mean mRNA copy number /i (see 
Eq. 20) may exist even in the absence of exotic DNA 



dynamics (jH (5] [10]) and could be attributed to the vari- 
ability in other extrinsic factors. We also predict that 
in contrast to proteins ([3]) the variation in intrinsic and 
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extrinsic factors both contribute significantly to noisy ex- 
pression of mRNA. Moreover, we directly probe the cell 
to cell variation in effective mRNA synthesis rate k and 
show that the coefficient of variation r/k plateaus at high 
expression levels (see bottom panel of Fig.|4|. 

We believe that the development presented here will 
serve as a framework to quantitatively parse the various 
contributors to noisy gene expression in more compli- 
cated situations e.g. regulated genetic circuits. 
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V. SUPPLEMENTARY MATERIALS 

A. Maximum entropy formalism 

The maximum entropy formalism allows one to esti- 
mate the probability distribution {pi\ of states {i} from 
limited information. Briefly, the maxEnt program for es- 
timating probabilities {pt} involves maximizing the en- 
tropy function S({pi}) subject to constraining the values 
of certain experimentally known observables (|23 | l24 p . For 
example, if (Xi), (X2), • ■ • , (Xn) are the mean values of 
observables X\, X%, . . . , Xn respectively, then the prob- 
abilities {pi} of states i are estimated by the maximizing 
the constrained objective function in Eq. [5j 

Here, {Xk} and a are Lagrange multipliers that en- 
sure that the experimental constraints are satisfied and 
that the probabilities are normalized. The entropy is a 
non-negative convex function of the probabilities and is 
usually defined as ([231 l2"o) 

Pi 

Here, p\ are the prior probabilities. 



B. Deriving the distribution of mRNAs 

Here's a detailed description of the maxEnt program 
to investigate the cell to cell variation in expression of 
mRNA molecules. 

As in the main text, consider the simplified reaction 
scheme 



a 



(5) 



DNA mRNA A 



(6) 

for transcription and degradation of mRNA molecules of 
a particular gene. 7 is the rate of transcription of mRNA 
molecules while 5 is the rate of degradation. The master 
equation describing the time evolution of the probability 
P(m|t, 7, S) of finding m mRNA molecules at time t is 
given by, 

dP(m\t) 



dt 



= 7PO - l\t) + (m + l)5P(m + l\t) 
- jP(m\t) - mSP(m\t). (7) 

For brevity, we have dropped the explicit dependence on 
7 and S. The steady state solution P(m|7, 8) of Eq.[7]with 
the initial condition that there were no mRNA molecules 
at t — is given by the Poisson distribution with param- 
eter k = j/5 (fT5]>. 



P(m|fc) 



(8) 



Since the parameter k varies from cell to cell, the distri- 
bution of mRNA copy numbers that is accessible exper- 
imentally is not P(m|fc) but in fact the marginal distri- 
bution, 



P(m) = ^P(m|fc) • P(k). 



(9) 



The maximum entropy formalism estimates P(k) from 
limited information. Let us denote the mean value of 
mRNA copy number at any particular value of k by 
(m)fc = fi. Experimentally, we know the population av- 
erage of mRNA copy number, 



P(k). 



(10) 



The maximum entropy estimate of P(k) is given by 
optimizing the constrained objective function, 

S [P(k)} +\iJ2 p {k) + ^J2( m )k ■ p (k) (11) 



Where, 



S = log 



m 
p°{k)' 



(12) 



Here P°(k) is the prior distribution, it is the probability 
P(fc) in the absence of any knowledge. It is customary 
to choose P°(k) as the entropic distribution 



where 



P°(k) oc exp(-X 3 S(k)) 



S(k) = -^P(m|/c)logP(m|fc). 



(13) 



(14) 



The entropic prior in Eq. 13 ensures that in the absence 
of any data, the less structured (broader) distribution 
functions P(m\k) are preferred over the peaked once. It 
can be shown that the choice of the entropic prior is 
equivalent to maximizing the entropy of the joint distri- 
bution P(m, k) rather than P(k). Interested readers are 
directed to ref. (JTHJ which provides a derivation for a 
general problem of the sort encountered here. 

Unfortunately, the entropy of the Poisson equation in 
Eq. [8] does not have a closed form but it grows logarith- 
mically for large k, 

S(k) « uiq +Wilogfe. (15) 

We use this form of the entropy in further calculations. 



Carrying out the optimization in Eq. 11 , we get the func- 
tional form, 



-i-*r [1 + /3] ' 



(16) 



The parameters a and /3 are related to the Lagrange mul- 
tipliers A2 and A3 respectively. 
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Consequently, the marginal distribution P(m) (See 
Eq|9} is given by the negative binomial distribution 
(Gamma distribution), 

a 1 +^(l + q)- 1 - m -^r[l + m + /3] 
P(m) = m!r[l + /3] ' (17) 

To simplify Eq. [16] and Eq. [TTJ let us calculate experi- 
mentally verifiable quantities from Eq. 17 (m) = /j,, the 
mean expression level of mRNA is given by, 



We thus parametrize Eq. 16 Eq. 17 in terms of [i and 
a and get 



P(fc; /i, a) 
P(m) 



e -ak Uafi—1 

a~^T [atfj] 

a a >"(l + a)- m - a ^T [m + a/j] 
m!T [afj] 



(21) 
(22) 



(i = (m) = 



1 + 



(18) 



The intrinsic noise 7/j, the extrinsic noise t^e, and the total 
noise ryx are given by, 



(19) 




(20) 



If we know P(m|fc) and P(k), the joint distribution 
P(m, fc) is given by 



P(m,/c) 



-(a+l)kum+an- 

a^^mlT [afj] 



(23) 



One can calculate an array of co-variances from 
the joint distribution. For example, the co-variance 

(mfc)-(m)(fc) • • , 

(m)W 18 s iven b ^ 



(mk) - (m) (k) _ J E m mfcP(m, k)dk - (/ kP{k)dk) (£ m mP(m)) _ 1 



i)(k) 



(JkP(k)dk) (£ m mP(m)) 



?7E 



(24) 



The co- variance can be experimentally estimated from 
a simultaneous measurement of the cell to cell vari- 
ation in the mRNA of a given gene and the ratio 
[RNAP]/[RNAase] of the copy numbers of RNA poly- 
merase holoenzyme and RNAase. 



C. Numerical simulations 

The synthesis and degradation of the mRNA of any 
given gene competes with the synthesis and degradation 
of all other co-expressed genes. Moreover, the cellular 
machinery that carries out these reactions itself com- 
prises of proteins and mRNAs and is subject to cell to 
cell variation. We devise a simple scheme to mimic the 
coupled dynamics of synthesis and degradation of the 
cellular machinery with the dynamics of synthesis and 
degradation of the mRNA of a given gene. 



I 

complex. Similarly, the rate of degradation S depends lin- 
early on the concentration [RNAase] of the proxy for the 
RNAase enzyme (7 = 7 [RNAP] and S = So [RNAase]). 
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RNAP 
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The transcription apparatus is represented by a single 
protein RNAP and the mRNA degradation apparatus is 
represented by a single protein RNAase. The rate of syn- 
thesis 7 of the given mRNA depends linearly on [RNAP] 
the concentration of the proxy for the RNA polymerase 



The dynamics of the synthesis and degradation of the 
mRNA of the given gene and RNAP and RNAase is prop- 
agated using the Gillespie's algorithm (j20|) for 2 • 10 8 
steps. Data is stored every 5000th step after an initial 
equilibration of 50000 steps. The initial concentrations 
of all species except the copy number of the each gene 



Parameters 

Parameter 


for the simulation 

Case 1 Case 2 Case 3 


7i 


9 n 


9 n 


9 n 


72 


9 n 




9 n 


70 


n q 

U . •) 


1 R 
1.1) 


1 n 


Kl 


n r i 


n r . 
u ■ o 


n r . 
u . o 


K2 


n r i 

u . o 


u . o 


fl r, 
u . u 


6l 


0.1 


0.1 


0.1 


52 


0.1 


0.1 


0.1 


So 


0.227 


0.5 


0.1 


A 1 


0.15 


0.65 


0.55 


A 2 


0.15 


0.65 


0.55 


[DNA] RNA p 


5 


5 


5 


[DNA] RNA aso 


5 


5 


5 


[DNAJggiic 


1 


1 


1 



TABLE I. The details of the parameters for the numerical simulation of mRNA synthesis. All rates are in s 1 and all copy 
numbers are integers. 



on the DNA at t = was set to 0. Table [T] gives the details of the conditions that were employed to construct 

the histograms (red points in Fig. 2 of the main text). 



